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Gerald  S.  Shedler 
IBM  Research  Laboratory 
San  Jose,  California 


ABSTRACT 

The  nonhomogeneous  Poisson  process  is 
a widely  used  model  for  a series  of  events 
(stochastic  point  process)  in  which  the 
"rate"  or  "intensity"  of  occurrence  of 
points  varies,  usually  with  time.  The 
process  has  the  characteristic  properties 
that  the  number  of  points  in  any  finite  set 
of  nonoverlapping  intervals  are  mutually 
independent  random  variables,  and  that  the 
number  of  points  in  any  of  these  intervals 
has  a Poisson  distribution.  In  this  paper 
we  first  discuss  several  general  methods 
for  simulation  of  the  one-dimensional  non- 
homogeneous Poisson  process;  these  include 
time-scale  transformation  of  a homogeneous 
(rate  one)  Poisson  process  via  the  inverse 
of  the  integrated  rate  function,  generation 
of  the  individual  intervals  between  points, 
and  generation  of  a Poisson  number  of  order 
statistics  from  a fixed  density  function. 

We  then  state  a particular  and  very 
efficient  method  for  simulation  of  nonhomo- 
geneous Poisson  processes  with  log-linear 
rate  function.  The  method  is  based  on  an 
identity  relating  the  nonhomogeneous  Poisson 
process  to  the  gap  statistics  from  a random 
number  of  exponential  random  variables  with 
suitably  chosen  parameters.  This  method 
can  also  be  used,  at  the  cost  of  program- 
ming complexity  and  some  memory,  as  the 
basis  for  a very  efficient  technique  for 
simulation  of  nonhomogeneous  Poisson  pro- 
cesses with  more  complicated  rate  functions 
such  as  a log-quadratic  rate  function. 


Finally,  we  describe  a simple  and  rela- 
tively efficient  new  method  for  simulation 
of  one-dimensional  and  two-dimensional  nci- 
homogeneous  Poisson  processes.  The  method 
is  applicable  for  any  given  rate  function 
and  is  based  on  controlled  deletion  of 
points  in  a Poisson  process  with  a rate 
function  that  dominates  the  given  rate  func- 
tion. In  its  simplest  implementation,  the 
method  obviates  the  need  for  numerical  in- 
tegration of  the  rate  function,  for  order- 
ing of  points,  and  for  generation  of  Poisson 
variates.  The  thinning  method  is  also 
applicable  to  the  generation  of  individual 
intervals  between  points,  as  is  required  in 
many  programs  for  discrete-event  simulations. 

1 . INTRODUCTION 

The  one-dimensional  nonhomogeneous 
(nonstationary)  Poisson  process  (see  e.g., 

[5,  pp.  28-29;  3,  pp.  94-101])  has  the 
characteristic  properties  that  the  numbers 
of  points  in  any  finite  set  of  nonoverlap- 
ping intervals  are  mutually  independent 
random  variables,  and  that  the  number  of 
points  in  any  interval  has  a Poisson  distri- 
bution. The  most  general  nonhomogeneous 
Poisson  process  can  be  defined  in  terms  of 
a monotone, nondecreas ing, right-continuous 
function  A(x)  which  is  bounded  in  any 
finite  interval.  Then  the  number  of  points 
in  any  finite  interval,  for  example  (0,Xg), 
has  a Poisson  distribution  with  parameter 
Uq  = A(Xg)  - A(0).  In  this  paper  we  assume 
that  A(x)  is  continuous,  but  not  necessarily 


1 


POISSON  PROCESSES  IN  NONSTATIONARY  SYSTEMS ...  Continued 


absolutely  continuous.  The  right  deriva- 
tive A(x)  of  A(x)  is  the  rate  function 
for  the  process;  A(x)  is  called  the  in- 
tegrated rate  function  and  has  the  inter- 
pretation that  for  x ^ 0,  A(x)  - A < 0 ) = 
E[N(x)],  where  N(x)  is  the  total  number 
of  points  in  (0,xj.  Note  that  A(x)  may 
jump  at  points  at  which  A(x)  is  not  abso- 
lutely continuous.  In  contrast  to  the 
homogeneous  Poisson  process,  i.e.,  A(x)  a 
constant  (usually  denoted  by  A) , the  in- 
tervals between  the  points  in  a one-dimen- 
sional nonhomogeneous  Poisson  process  are 
neither  independent  nor  identically  dis- 
tributed. 

Applications  of  the  one-dimensional 
nonhomogeneous  Poisson  process  include 
modelling  of  the  incidence  of  coal-mining 
disasters  [5],  the  arrivals  at  an  inten- 
sive care  unit  [12],  transaction  processing 
in  a data  base  management  system  [15],  oc- 
currences of  major  freezes  in  Lake  Constance 
[23],  and  geomagnetic  reversal  data  [22]. 
The  statistical  analysis  of  trends  in  a 
one-dir.^nsional  nonhomogeneous  Poisson  pro- 
cess, based  on  the  assumption  of  an  expo- 
nential polynomial  rate  function,  is  dis- 
cussed by  [4,  5,  12  and  15]. 

One-dimensional  nonhomogeneous  Poisson 
processes  are  often  used  as  models  for 
event  streams  when  there  is  gross  inhomo- 
geneity in  a system,  e.g.,  time  of  day 
effect  or  long-term  growth  in  use  of  a 
facility.  It  is  important  to  be  able  to 
simulate  these  processes  since  analytic 
results  are  difficult  to  obtain.  This  is 
particularly  true  in  the  context  of  queue- 
ing systems;  see  e.g.,  [19].  The  methods 
given  here  for  simulation  of  the  one-dimen- 
sional nonhomogeneous  Poisson  process  have 
application,  for  example,  to  study  of  the 
length  of  a queue  at  a toll  booth  at  a 
time  corresponding  to  the  peak  traffic 
time,  or  to  study  of  the  arrivals  at  an  in- 
tensive care  unit  where  the  probability  of  a 


bed  being  free  at  a time  corresponding  to  the 
peak  of  arrivals  from  afternoon  operations 
is  of  interest.  Note  that  in  simulations 
of  nonhomogeneous  systems  of  this  kind, 
estimates  of  measures  of  system  behavior 
will  be  based  on  multiple  replications. 

The  two-dimensional  homogeneous  Poisson 
process  (of  rate  A > 0)  is  defined  by  the 
properties  that  the  numbers  of  points  in 
any  finite  set  of  nonoverlapping  regions 
having  areas  in  the  usual  geometric  sense 
are  mutually  independent,  and  that  the  num- 
ber of  points  in  any  region  of  area  A has 
a Poisson  distribution  with  mean  *A;  see 
e.g.,  [10,  pp.  31-32],  Note  that  the  num- 

ber of  points  in  a region  R depends  on  its 
area,  but  not  on  its  shape.  The  homogeneous 
Poisson  process  arises  as  a limiting  two- 
dimensional  point  process  with  respect  to 
a number  of  limiting  operations;  cf . , [7,8]. 

Properties  of  the  process  are  given  in  [18]. 
Applications  of  the  two-dimensional  homo- 
geneous Poisson  process  to  problems  in 
ecology  and  forestry  have  been  discussed  in 
[24]  and  [9].  The  model  also  arises  in 
connection  with  naval  search  and  detection 
problems . 

The  two-dimensional  nonhomogeneous 
Poisson  process  is  characterized  by  a (con- 
tinuous) positive  rate  function  (x,y). 
Applications  of  the  two-dimensional  nonhomo- 
geneous Poisson  process  include  problems  in 
forestry  as  well  as  naval  search  and  detec- 
tion. The  detection  and  statistical  analy- 
sis of  trends  in  the  two-dimensional  non- 
homogeneous Poisson  process  is  discussed  in 
[21]  . 

2.  SIMULATION  OF  THE  ONE-DIMENSIONAL 
NONHOMOGENEOUS  POISSON  PROCESS 

There  are  a number  of  methods  for  sim- 
ulating one-dimensional  nonhomogeneous 
Poisson  process  which  we  review  briefly. 
Time-scale  transformation  of  a homogeneous 
(rate  one)  Poisson  process  via  the  inverse 
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of  the  integrated  rate  function  A(x)  con- 
stitutes a first  general  method;  cf.  [3, 
pp.  96-97],  This  method  is  based  on  the 
result  that  , X2 , ...  are  the  points 
in  a nonhomogeneous  Poisson  process  with 
continuous  integrated  rate  functions  A(x) 
if  and  only  if  X|  = AfX^,  X£  = A(X2),..., 
are  the  points  in  a homogeneous  Poisson 
process  of  rate  one.  The  time-scale  trans- 
formation method  is  a direct  analogue  of 
the  inverse  probability  integral  transfor- 
mation method  for  generating  (continuous) 
nonuniform  random  numbers.  For  many  rate 
functions,  inversion  of  A(x)  is  not 
simple  and  must  be  done  numerically;  cf . , 
[6]  and  [20].  The  resulting  algorithm  for 
simulation  of  the  nonhomogeneous  Poisson 
process  may  be  far  less  efficient  than 
simulation  based  on  other  methods. 

A second  general  method  for  simulating 
a one-dimensional  nonhomogeneous  Poisson 
process  with  integrated  rate  function  A(x) 
is  to  generate  the  intervals  between  points 
individually,  an  approach  which  may  seem 
rrore  natural  in  the  event  scheduling  ap- 
proach to  simulation.  Thus,  given  the 
points  X1  = x^,  X2  = x2 , . . . , X^  = x^ , -fith 
Xi  < X2  < •••  < Xi»  the  interval  to  the 
next  point,  X.+1  - X.,  is  independent  of 

xi_i  and  1133  distribution  function 
F(x)  = 1 - exp[-{A(xi  + x)  - A(x^)}].  It 
is  possible  to  find  the  inverse  distribu- 
tion function  F—  ■*"(•),  usually  numerically, 
and  generate  Xi+1  - according  to 

Xi+1  ” xi  = where  IT  is  a uni- 

form random  number  on  the  interval  (0,1). 
Note,  however  that  this  not  only  involves 
computing  the  inverse  distribution  function 
for  each  interval  X^+1  - x^»  but  that  each 
distribution  has  different  parameters  and 
possibly  a different  form.  An  additional 
complication  is  that  X^+^  - x^  is  not 
necessarily  a proper  random  variable,  i.e. 
there  may  be  positive  probability  that 
X^+j  - X^^  is  infinite.  It  is  necessary 
to  take  this  into  account  for  each  interval 
X^+^  - X^  before  the  inverse  probability 


integral  transformation  is  applied.  The 
method  is  therefore  very  inefficient  with 
respect  to  speed,  more  so  than  the  time- 
scale  transformation  method. 

In  a third  method,  simulation  of  a 

non-homogeneous  Poisson  process  in  a fixed 

interval  (0,  xQ ] can  be  reduced  to  the 

generation  of  a Poisscn  number  of  order 

statistics  from  a fixed  density  function 

by  the  following  result;  cf. , [5,  p.  45], 

If  X^ , X2,  ...  , Xn  are  the  points  of  the 

nonhomogeneous  Poisson  process  in  (0,  xQ ] , 

and  if  N(Xq)  = n,  then  conditional  on 

having  observed  n ( > 0)  points  in  (0,xQ), 

the  X.  are  distributed  as  the  order  sta- 
l 

tistics  from  a sample  of  size  n from  the  dis- 
tribution function {A(x)-A(0)}/iA(x-)-A(0) ; , 
defined  for  0 < x <_  x^ . Simulation  of  the 
nonhomogeneous  Poisson  process  based  on 
order  statistics  is  in  general  more  effi- 
cient (with  respect  to  speed)  than  either 
of  the  previous  two  methods.  Of  course,  a 
price  is  paid  for  this  greater  efficiency. 
First,  it  is  necessary  to  be  able  to  gen- 
erate Poisson  variates,  and  second,  more 
memory  is  needed  than  in  the  interval-by- 
interval method  to  store  the  sequence  of 
points.  Enough  memory  must  be  provided  so 
that  with  very  high  probability  the  random 
number  of  points  generated  in  the  interval 
can  be  stored.  Recall  that  the  number  of 
points  in  the  interval  (0,Xq]  has  a Poisson 

distribution  with  mean  u0  = .'.(x0)  - A ( 0 ) . 

1/2 

Memory  of  size,  e.g.,  Mq + 4uq  will  en- 
sure that  overflow  will  occur  on  the  aver- 
age in  only  1 out  of  approximately  40,000 
realizations.  This  probability  is  small 
enough  so  that  in  case  of  overflow,  the 
realization  of  the  process  generally  can 
be  discarded. 

We  now  summarize  several  recently 
developed  methods  for  simulating  one-  and 
two-dimensional  nonhomogeneous  Poisson  pro- 
cesses. These  methods  are  discussed  in 
greater  detail  in  [14,  16,  17]. 
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3.  SIMULATION  USING  GAP  STATISTICS 

In  a previous  paper  [14],  we  have  con- 
sidered the  simulation  of  nonhomgeneous 
Poisson  processes  with  degree-one  exponen- 
tial polynomial  rate  function,  i.e., 

A(x)  of  the  form 

\ (x)  = exp{yQ  + y^x}  = X exp{y^x],  y^O.  (1) 

The  rate  function  (1)  is  the  simplest  of 
a general  family  of  log-linear  rate  func- 
tions, i.e.,  rate  functions  whose  loga- 
rithms are  linear  in  the  coefficients  [4; 

12]  which  are  useful  in  analyzing  nonhomo- 
generous  Poisson  processes.  The  rate 
function  (1)  represents  a situation  in 
which  the  rate  is  monotonically  increasing 
or  decreasing  depending  on  whether  y^  is 
greater  than  or  less  than  zero,  with  y^ 
equal  to  zero  giving  a homogeneous  Poisson 
process.  The  case  where  y^  is  less  than 
zero  and  the  case  where  y^  is  greater 
than  zero  are  quite  distinct;  in  the  first 
situation  X(x)  •*  0 as  x -*■  ",  and  in  the 
second,  A(x)  -*•  " as  x ».  Moreover, 
when  y ^ is  less  than  0,  the  intervals 
between  events  are  not  proper  random  vari- 
ables since  there  is  a nonzero  probability 
that  there  is  no  event  after  any  fixed 
point  x. 

In  [14]  a method  for  simulating  the 
nonhomogeneous  Poisson  process  is  given 
based  on  an  identity  relating  the  nonhomo- 
geneous Poisson  process  with  rate  function 
(1)  to  the  gap  statistics  from  a random 
number  of  exponential  random  variables  with 
suitably  chosen  parameters.  This  method 
avoids  costly  ordering  and  taking  of  loga- 
rithms required  by  direct  simulation  methods 
and  is  more  efficient  than  time-scale  trans- 
formation of  a homogeneous  Poisson  process 
via  the  inverse  of  the  integrated  rate 
function  A(x). 

Simulation  of  the  one-dimensional 
nonhomogeneous  Poisson  process  in  a 
fixed  interval  (0,xn]  is  more 
natural  than  simulation  for  a fixed 
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number  of  events  since  time,  not  serial 
number,  is  the  basic  parameter  of  the  in- 
homogeneity discussed  here.  Thus  the  gap 
statistics  algorithm  for  simulation  of  the 
nonhomogeneous  Poisson  process  generates 
the  sequence  of  times-to-events  in  a fixed 
interval.  Although  such  a method  requires 
more  memory  than  successive  generation  of 
individual  times  until  the  next  event,  it 
is  far  more  efficient. 

We  now  state  the  method  of  Lewis  and 
Shedler  [14]  for  simulation  via  gap  statis- 
tics of  the  one-dimensional  nonhomogeneous 
Poisson  process  with  rate  function  (1) . 

This  scheme,  which  is  particular  to  the 
degree-one  exponential  polynomial  rate 
function,  can  use  standard  packages  for  ex- 
ponential random  numbers  (e.g.,  [13])  and 
obviates  the  need  tor  ordering  of  the  random 
numbers.  It  is  based  on  the  result  (see 
[25])  that  the  gap  process  associated  with 
a Poisson  distributed  (parameter  ->/y^  >0) 
number  of  exponential  (parameter  8 = -y^>0) 
gap  statistics  is  a nonhomogeneous  Poisson 
process  with  rate  function  A(x)  = 1 expty^x) 
on  (0,~).  Efficient  methods  for  genera- 
tion of  Poisson  random  numbers  for  which 
the  generation  time  does  not  increase  pro- 
portionally with  the  mean  are  given  by  [1,2] 
and  [11]. 

Assuming  the  availability  of  a source 
of  unit  exponential  random  numbers  E^,E,, 
...,  obtained  by  logarithms  or  by  other 
methods,  the  resulting  algorithm  for  gener- 
ating the  events  in  the  nonhomogeneous 
Poisson  process  is  as  follows. 

Algorithm  1.  Gap  Statistics  Technique 

(YX  < 0). 

1.  Generate  m as  a Poisson  random  number 
with  parameter  -X/y..  If  m = 0,  exit; 
there  are  no  events  in  (0,XgJ. 

2.  For  m > 0,  if  E^/fBm)  is  greater 
than  Xq , exit;  there  are  no  events  in 
(0,Xg).  Otherwise,  set  E^/(6m)  equal 
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extended  to  more  complex  rate  functions 
than  the  degree-two  exponential  polynomial. 


3.  If  E2/{0  (m-1)  } + > xQ,  then  return 

and  exit.  Otherwise,  set  it  equal 

to  X2. 

4.  Continue,  possibly  for  m times.  If 

V6  + Vl  > x0'  return  X1'X2 Xm-1 

and  exit.  Otherwise,  set  this  equal 

to  X , return  X.,X~, ,X„  and  exit. 

m i ^ m 

The  case  y1  > 0 is  handled  in  the 
same  way  as  < 0 by  using  a time-rever- 

sal technique,  as  follows.  Simulate  ac- 
cording  to  Algorithm  1 with  A(x)  = A 
x exp{y  *x]  , where  A*  = exp{yQ  + y^g} 
and  y*  = - Yj^-  The  output  of  Algorithm  1 
is  a sequence  X* , X2 , . . . ,Xn . Then  set 

X1  = X0  " Xn ' X2  = x0~  Xn-1  - • • • »xn  = x0  ~ x*  • 

These  X.  are  the  required  events  in  the 

nonhomogeneous  Poisson  process  for  y^  > 0. 

Lewis  and  Shedler  [16]  consider  the 
simulation  in  a fixed  interval  of  the  non- 
homogeneous Poisson  process  with  degree- 
two  exponential  polynomial  rate  function 

A (x)  = exp{aQ  + a^x + a2x2} , a2 ? 0 (2) 

the  case  a2  = 0 giving  the  degree-one 
exponential  polynomial  rate  function. 

Again,  the  case  where  a2  is  less  than 
zero  and  the  case  where  a2  is  greater 
than  zero  are  distinct.  The  simulation 
method  give/,  is  based  on  representation  of 
the  process  as  a superposition  of  two  inde- 
pendent nonhomogeneous  Poisson  processes, 
one  of  which  has  a fitted  rate  function  of 
the  form  (1) ; simulation  of  the  latter 
process  is  accomplished  via  the  gap  statis- 
tics algorithm  [14].  A rejection-accep- 
tance technique  is  used  to  generate  the 
other,  more  complex,  nonhomogeneous  Poisson 
process.  The  resulting  algorithm  is  more 
efficient  than  time-scale  transformation  of 
a homogeneous  Poisson  process;  see  [20]. 
This  method  can  be  improved  by  using  the 
thinning  algorithm  given  in  the  next  sec- 
tion to  simulate  the  second  nonhomogeneous 
Poisson  process.  The  method  can  also  be 


4.  SIMULATION  OF  NONHOMOGENEOUS 
POISSON  PROCESSES  BY  THINNING 

In  this  section  we  describe  a new 
method  [17]  for  simulating  a nonhomogeneous 
Poisson  process.  The  method  is  not  only 
conceptually  simple,  but  is  also  computa- 
tionally simple  and  relatively  efficient. 

In  fact,  at  the  cost  of  some  efficiency, 
the  method  can  be  applied  to  simulate  the 
given  nonhomogeneous  Poisson  process  with- 
out the  n^ed  for  numerical  integration  or 
routines  for  generating  Poisson  variates ■ 
Used  in  conjunction  with  the  special  methods 
given  by  Lewis  and  Shedler  [14,16],  the 
method  can  be  used  to  simulate  quite  ef- 
ficiently nonhomogeneous  Poisson  processes 
with  rather  complex  rate  function,  in 
particular,  combinations  of  long-term  trends 
and  fixed-cycle  effects.  The  method  is 
also  easily  extended  to  the  problem  of  simu- 
lating the  two-dimensional  nonhomogeneous 
Poisson  process  (see  Section  5) , and  of 
simulating  conditional  and  doubly  stochastic 
Poisson  processes. 


Simulation  of  a nonhomogeneous  Poisson 
process  with  general  rate  function  *(x) 
in  a fixed  interval  (0,Xg]  can  be  based 
on  thinning  of  a nonhomogeneous  Poisson 

1c 

process  with  rate  function  A (x)  > (x)  . 

The  main  result  [17,  Theorem  1]  is  that  if 

* * * 

X^ ,X2 , . . . , j are  the  points  of  the 
0 * 

process  with  rate  function  A (x)  in  the 

* 

interval  (0,Xg]  and  if  the  point  XV  is 
deleted  with  (independent)  probability 
1 - X(X*)/A*(X*)  , then  the  remaining  points 
form  a nonhomogeneous  Poisson  process  with 
rate  function  A(x)  in  the  interval  (0,Xq]. 

This  result  is  the  basis  for  the  method 
of  simulating  one-dimensional  nonhomogeneous 


Poisson  processes  on  an  interval  (0,xQ] 
given  by  Algorithm  2. 


5 


POISSON  PROCESSES  IN  NONSTATIONARY  SYSTEMS ...  Continued 


Algorithm  2.  One-Dimensional  Nonhomogen- 

eous  Poisson  Process  (Thinning) 

1.  Generate  points  in  the  nonhomogeneous 

* 

Poisson  process  {N  (x) } with  rate 

* 

function  X (x)  in  the  fixed  interval 

(0,x_].  If  the  number  of  points  gen- 
0 * * 
erated,  n , is  such  that  n =0,  exit; 

there  are  no  points  in  the  process 

{N (x) } . 

* * 

2.  Denote  the  (ordered)  points  by  X.,X.., 

* ■*.  * 
...»  Xn*.  Set  i = 1 and  k = 0 . 

3.  Generate  IK,  uniformly  distributed 
between  0 and  1.  If  IK£A(X^)/X  (X^)  , 


set 

k 

equal  to 

k+1 

and 

*k 

= X*. 
1 

4. 

Get 

i 

equal  to 

i+1 

. If 

i < 

* 

. n ' 

go 

to  3. 

5. 

Return 

X.  , X, , . • . 

'Xn' 

where 

n 

= x, 

and 

also  n. 

The  method  of  thinning  of  Algorithm  2 
is  essentially  the  obverse  of  the  condi- 
tional method  of  Section  2 using  condition- 
ing and  acceptance-rejection  techniques  to 
generate  the  random  variables  with  density 
function  X(x)/{A(x)  - A ( 0 ) > [16,  Algorithm 
3].  The  differences  are  subtle,  but  compu- 
tationally important.  In  the  acceptance- 
rejection  method,  it  is  first  necessary  to 
generate  a Poisson  variate  with  mean 
= A(xq)  - A(0),  and  this  involves  an 
integration  of  the  rate  function  X(x). 

Then  the  Poisson ( Ug)  number,  n,  of  vari- 
ates generated  by  acceptance-rejection  must 
be  ordered  to  give  xi  ,X2 ' ‘ ‘ ' Xn' 

In  the  simplest  form  of  the  method  of 
thinning,  X (x)  is  taken  equal  to 
X * maxo  < x < x ' so  that»  for  instance 

• * • 

the  points  xi »X2 ' • • • ,Xn*  can  be  generated 
by  cumulating  exponential ( A * ) variates  until 
the  sum  is  greater  than  xQ  (cf.  [14,  Al- 
gorithm 1].  The  generated  points  are  then 
thinned.  No  ordering,  no  integration  of 
A(x)  and  no  generator  of  Poisson  variates 


is  required.  Of  course  for  both  algorithms 
to  be  efficient,  computation  of  X(x)  and 

h 

X (x)  must  be  easy  relative  to  computation 
of  the  inverse  of  A(x). 

For  the  thinning  algorithm  (as  well  as 
the  algorithm  based  on  conditioning  and 
acceptance-rejection)  efficiency  as  measured 
by  the  number  of  points  deleted  is  propor- 
tional to 

UQ/U*Q  = f A (Xq ) - A (0)  }/{ A*  (xQ ) - A*  (0)  ) ; 


this  is  the  ratio  of  the  areas  between  0 

* 

and  xn  under  X(x)  and  X (x) . Thus, 
X (x)  should  be  as  close  as  possible  to 
X(x)  consistent  with  ease  of  simulating 
the  nonhomogeneous  Poisson  process 
{N* (x)  :x  > 0}  . 


It  is  important  to  note  that  the  method 
of  thinning  can  be  used  to  generate  individ- 
ual intervals  between  events  occurring  in 
(0,xQ]  if  > (x ) is  bounded  on  (0,xQl. 

The  resulting  algorithm  is  not  only  useful 
in  the  event  scheduling  approach  to  simula- 
tion, but  also  in  the  generation  of  con- 
ditional Poisson  processes.  Informally, 

the  ’bne  at  a time"  thinning  algorithm  is 

* * 

as  follows.  If  X (x)  = A ^maxo  < x < x X ' 

then  the  ith  interval  X.  - X.  , is  ob- 

l i-l 

tained  by  generating  and  cumulating  expo- 

* ★ 

nentialfX*)  random  numbers  E,  ,,  E.  ,,... 
until,  for  the  first  time, 


Ui.j  i X(Xi-l  + Ei,l  +-*'+  Ei,j)/X 


i — 1,2,...;  j — 1,2,..., 


where  the  U.  . are  independent,  uniform 
t , 3 

(0,1)  random  numbers. 


5.  SIMULATION  OF  TWO-DIMENSIONAL 
HOMOGENEOUS  POISSON  PROCESSES 

Recall  that  the  two-dimensional  homo- 
geneous Poisson  process  (of  rate  X > 0)  nas 
the  characteristic  properties  that  the  num- 
bers of  points  in  any  finite  set  of 
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nonoverlapping  regions  having  areas  in  the 
usual  geometric  sense  are  mutually  inde- 
pendent, and  that  the  number  of  points  in 
any  region  of  area  A has  a Poisson  dis- 
tribution with  mean  AA. 

In  considering  the  two-dimensional 
homogeneous  Poisson  process,  projection 
properties  of  the  process  depend  quite 
critically  on  the  geometry  of  tne  regions 
considered.  These  projection  properties 
are  simple  for  rectangular  and  circular 
regions,  and  make  simulation  of  the  homo- 
geneous process  quite  easy.  We  consider 
here  the  case  of  a rectangular  region.  The 
following  result  forms  the  basis  for  simu- 
lation of  the  two-dimensional  homogeneous 
Poisson  process  of  rate  X in  a fixed 
rectangle  R = {(x,y):0£x£x0,0<y£  yQ  }. 

If  (X^Y^,  (X2,Y2)  , . . . , (Xjj.Yjj)  denote 

the  positions  of  the  points  of  the  process 
in  R,  labelled  so  that  X^  < X2  < " * , then 
X^,X2,...,Xn  form  a one-dimensional  homo- 
geneous Poisson  process  on  0 £ x £ xQ  of 
rate  Ay^;  if  tne  points  are  relabelled 
(X^,Y^),  (X2, , Y2 ) , ...  , (XjJ. , YjJj)  so  that 

Y[  < Y2  <•••<  Y^,  then  Y|,Y^,...,Y^  form 

a one-dimensional  homogeneous  Poisson  pro- 
cess on  0 < Y £ Yq  of  rate  Axg. 

We  state  next  conditional  properties 
of  the  Poisson  process  in  a rectangle.  The 
important  thing  to  note  is  that  although 
the  processes  obtained  by  projection  of  the 
points  onto  the  x and  y axes  are  not 
independent,  there  is  a type  of  conditional 
independence  which  makes  it  easy  to  simu- 
late the  two-dimensional  process.  Thus, 
conditional  on  having  observed  n > 0 points 
(X1,Yi),  (X2,Y2) .... , (Xn,Yn)  in  R,  labelled 
so  that  Xx  < X2  <•••<  x , the  X1,X2,...,Xn 
are  uniform  order  statistics  on  0 £ x £ xQ , 
and  the  Y^,Y2,...,Y  are  independent  and 
uniformly  distributed  on  0 < y £ yQ , in- 
dependent of  the  X^ . 

From  these  two  results,  the  following 
simulation  procedure  is  obtained. 


Algorithm  3.  Two-Dimensional  Homogeneous 

Poisson  Process  in  a Rectangle 

1.  Generate  points  in  the  one-dimensional 
homogeneous  Poisson  process  of  rate 
Ay^  on  (0,Xg],  If  the  number  of 
points  generated,  n,  is  such  that  n=0, 
exit;  there  are  no  points  in  the 
rectangle . 

2.  Denote  the  points  generated  by 

12  n 

3.  Generate  Y^,Y2,...,Y  independent, 
uniformly  distributed  random  numbers 
on  (0 ,yQ  ] • 

4.  Return  (Xj^Y^  , (X2  ,Y2)  , . . . , (Xn  ,Yn>  as 
the  coordinates  of  the  two-dimensional 
homogeneous  Poisson  process  in  the 
rectangle,  and  n. 

Note  that  generation  of  the  points  X^,X2, 

. ..,X  in  steps  1 and  2 can  be  accomplished 
by  cumulating  exponential ( Ay  ) random 
numbers.  Alternatively,  after  generating  a 
Poisson  random  number  N = n (with  parameter 
Ax0y0),  n independent,  uniformly  distributed 
random  numbers  on  (0,Xg]  can  be  ordered; 
see  [14,  p.  502}. 

The  basis  for  another  algorithm  for  sim- 
ulation of  the  two-dirensional  homogeneous 
Poisson  process  in  a rectangle  is  the  follow- 
ing corollary.  Denote  the  Poisson  points  by 
(Xf, Yf) , (X2 , Y2> , . . . , where  the  index  dees 
not  necessarily  denote  an  ordering  on  either 
axis.  Conditionally,  the  pairs  (X^.Y^), 

(X2 , Y2) , . . . , (X^ , Y ) are  independent  random 
variables,  and  furthermore,  for  each'  pair 

(X.,Y.),  X.  is  distributed  uniformly  between 
ill 

0 and  Xg , independently  of  Y^ , which  is 
uniformly  distributed  between  0 and  y^. 

Direct  generation  of  homogeneous  Poisson 
points  in  noncircular  or  nonrectangul ar 
regions  is  difficult.  The  processes  ob- 
tained by  projection  of  the  points  on  the 
two  axes  are  non homogeneous  Poisson  processes 
with  complex  rate  functions  determined  by 
the  geometry  of  the  region.  However,  the 
conditional  independence  which  is  found  in 
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circular  and  rectangular  regions  for  the 
processes  on  the  two  axes  is  not  present. 

In  particular,  given  that  there  are  n 
points  ^xi»Yj_)  »•••.»  (X^Y  ) in  a nonrec- 
tangular  region,  the  pairs  (X^,Y^)  are 
mutually  independent,  but  is  in  gen- 

eral not  independent  of  Y.,  i = l,2,...,n. 
Therefore  it  is  simpler  to  enclose  the 
region  in  either  a circle  or  a rectangle, 
generate  a homogeneous  Poisson  process  in 
the  enlarged  area  , and  subsequently  ex- 
clude points  outside  of  the  given  region. 

6 . SIMULATION  OF  TWO-DIMENSIONAL 

NONHOMOGENEOUS  POISSON  PROCESSES 

The  two-dimensional  nonhomogeneous 
Poisson  process  {N(x,y):x£0,  y£0}  is 
specified  by  a positive  rate  function 
X(x,y)  which  for  simplicity  is  assumed 
here  to  be  continuous.  Then  the  process 
{N(x,y)}  has  the  characteristic  properties 
that  the  numbers  of  points  in  any  finite 
set  of  nonoverlapping  regions  having  areas 
in  the  usual  geometric  sense  are  mutually 
independent,  and  that  the  number  of  points 
in  any  such  region  R has  a Poisson  dis- 
tribution with  mean  A(R);  here  A ( R)  de- 
notes the  integral  of  X(x,y)  over  R, 

i.e. , over  the  entire  area  of  R. 

The  basic  result  of  Section  4 on 
thinning  of  one-dimensional  nonhomogeneous 
Poisson  processes  generalizes  to  two-dimen- 
sional nonhomogeneous  Poisson  processes. 

* 

Thus,  suppose  that  X(x,y)  <_  X (x,y)  in  a 
fixed  rectangular  region  of  the  plane.  If 

a nonhomogeneous  Poisson  process  with  rate 

* 

function  X (x,y)  is  thinned  according  to 

* 

Kx,y)/X  (x,y)  (i.e.,  each  point  (X^Y^ 

is  deleted  independently  if  a uniform  (0,1) 
random  number  U.  is  greater  than 

t 

X(X^,Y^)/X  (X£*Y^)),  the  result  is  nonhomo- 
geneous Poisson  process  with  rate  function 
A(x,y) . 

The  nonhomogeneous  Poisson  process  with 
rate  function  X(x,y)  in  an  arbitrary  but 


fixed  region  R can  be  generated  by  en- 
closing the  region  R either  in  a circle  [17], 
or  a rectangle,  and  applying  Alaorithm  3.  The 

following  procedure  assumes  that  the  region 

★ 

R has  been  enclosed  in  a rectangle  R , 

* 

and  that  X = max( X (x,y) :x,  y t R)  has 

been  determined;  here  the  bounding  orocess 

★ 

is  homogeneous  with  rate  X in  the  rectan- 
gle R . 

Algorithm  4.  Two-Dimensional  Nonhomoaeneous 
Poisson  Process  (Thinning) 

1.  Using  Algorithm  2,  generate  points  in 

the  homogeneous  Poisson  orocess  of  rate 
* * 

X in  the  rectangle  R . If  the  num- 

★ ★ 

ber  of  points,  n , is  such  that  n =0, 
exit;  there  are  no  points  in  the  non- 
homogeneous Poisson  process. 

* 

2.  From  the  n points  generated  in  1, 
delete  the  points  that  are  net  in  R, 
and  denote  the  remaining  points  by 

(xJ-yJ).  U2S2) <VV  with 

X*  < X*  <•••<  X . Set  i = 1 and  k = 0. 
12  m 

3.  Generate  U^  uniformly  distributed 
between  0 and  1.  If  U^  £X(X^,Y^)/>  , 

* * 

set  k = k+1,  X^  = and  Y^  = Y^. 

* 

4.  Set  i equal  to  i+1.  If  i £ m , go 
to  3 . 

5.  Return  (X^Yj)  , (Xj.Yj),...,  (Xn,Yn>, 

where  n = k and  n. 

7.  CONCLUSION 

We  have  summarized  previously  known 
general  methods  for  simulating  nonhomogen- 
eous Poisson  processes  in  one  dimension. 

In  addition,  we  have  described  the  sirrrole 
and  efficient  new  methods  of  Lewis  and 
Shedler  for  simulating  nonhomogeneous 
Poisson  processes  in  one  and  two  dimensions. 
Extensions  of  the  thinning  algorithm  to  the 
simulation  of  homogeneous  or  nonhomogeneous 
conditional  or  doubly  stochastic  Poisson 
processes  will  be  described  elsewhere. 
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